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The computational study of conformational transitions in RNA and proteins with atomistic molec¬ 
ular dynamics often requires suitable enhanced sampling techniques. We here introduce a novel 
method where concurrent metadynamics are integrated in a Hamiltonian replica-exchange scheme. 
The ladder of replicas is built with different strength of the bias potential exploiting the tunability of 
well-tempered metadynamics. Using this method, free-energy barriers of individual collective vari¬ 
ables are significantly reduced compared with simple force-field scaling. The introduced methodology 
is flexible and allows adaptive bias potentials to be self-consistently constructed for a large number 
of simple collective variables, such as distances and dihedral angles. The method is tested on alanine 
dipeptide and applied to the difficult problem of conformational sampling in a tetranucleotide. 


Introduction 

Biomolecular dynamics involves a wide range of timescales ranging from bond fluctuations to slow large-scale 
motions. m Molecular dynamics (MD) with accurate force fields can in principle be used as a virtual microscope to 
investigate these motions at atomistic resolution. [1] However, its applicability to problems such as folding or conforma¬ 
tional transitions in proteins and RNA is limited by the fact that only short time scales /us) are directly accessible by 
straightforward simulation. In spite of the development of ad hoc hardware, [5] many relevant processes are still out of 
reach for accurate atomistic modeling. Several different techniques have been developed in the last decades to address 
this issue. These techniques can be roughly classified in two groups. [H] In the first group, inspired by annealing, [7] 
ergodicity is achieved by increasing the temperature [8l [9] or by artificially modifying the Hamiltonian. 1 1 OlfT^ In 
methods of this group a ladder of replicas with different degree of ergodicity is often employed. Swaps of coordinates 
between neighboring replicas are periodically attempted and accepted or rejected with a Metropolis criterion. These 
methods are usually referred to as replica exchange molecular dynamics (REMD). Whereas temperature is certainly 
the most adopted control parameter, temperature REMD (T-REMD) is computationally demanding for solvated 
systems, because replica spacing and interchange probability depends on the system size.[T3] Moreover, T-REMD is 
ineffective on entropic barriers. [HHS] Scaling of portions of the Hamiltonian (H-REMD) is a common alternative 
and could have a better convergence behavior for large systems. T-REMD and H-REMD can also be combined, by 
integrating both schemes on each renlicapT] or in a multidimensional framework. |I8j The second group of enhanced 
sampling techniques includes methods based on importance sampling, where suitable collective variables (GVs) are a 
priori selected and biased. This class has its root in umbrella sampling method, [TS] and includes local elevation. ^0] 
conformational hooding, [5T] adaptive biasing force (ABF), [52] and metadynamics. [23 I21| These methods are very 
efficient but require large a priori information. With the notable exception of bias-exchange metadynamics, |25j 
approaches based on importance sampling have been traditionally applied to a small number of CVs at a time, due 
to the difficulties in building a history dependent potential in a high-dimensional CV space. For many systems it is 
difficult to find a small number of effective CVs that describe the slow degrees of freedom and one often has to resort 
to expensive methods of the first class. For instance, conformational transitions in unstructured oligonucleotides have 
been studied with different REMD schemes. [MUSI HZ] In these studies, the generation of a converged conformational 
ensemble was proven a cumbersome task, even when a high number of replicas and tens of pLS of simulated time were 
employed. 

Methods bridging between these two classes can be designed by biasing a large number of local CVs (e.g. dihedral 
angles), so as to avoid the complication of designing ad hoc CVs. For example, Straatsmaa and McCammon[55] 
introduced a technique where bias potentials acting on dihedrals was used in a simulated annealing protocol. In 
that work a bias was fitted to the potential of mean force of backbone dihedrals and then used to quickly optimize 
the structure of a polypeptide. Zacharias and collaborators [23 HO] used a similar technique to build bias potentials 
for dihedrals to be employed in H-REMD and successfully applied it to study conformational changes in proteins. 
However, these potentials were fitted on a model system and cannot account for the specific identity of each residue 
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and for the cross-talk between correlated dihedrals. For nucleic acids, where the complex backbone does not allow 
a straightforward application of this technique, penalty potentials centered on the stable rotamers were manually 
selected with a procedure that seems difficult to generalize. [?TH55] 

In this paper we propose to use concurrent well-tempered metadynamics (53] (WT-MetaD) to build bias potentials 
acting on a large number of local CVs. We then show how to integrate this approach in a H-REMD scheme, exploiting 
the replica ladder to obtain unbiased conformations. In WT-MetaD the compensation of the underlying free-energy 
landscape is modulated by the so-called bias factor 7 . We here change this parameter across the replica ladder, 
adjusting the ergodicity of each replica. The final bias could be also used as a static potential so as to completely 
eliminate any non-equilibrium effect. Since the effect of the bias is that of keeping the chosen CVs at an effectively 
higher temperature, we refer to the introduced method as replica exchange with collective-variable tempering (RECT). 
The method is first tested on alanine dipeptide in water and then applied to the conformational sampling of a RNA 
tetranucleotide where it outperforms dihedral-scaling REMD and plain MD. The chosen tetranucleotide is a very 
challenging system that has been recently studied with different variants of REMD. 


Methods 

In this section we show how to use WT-MetaD as an effective method to build bias potentials that allow barriers to be 
easily crossed. One of the input parameters of well-tempered metadynamics is a boosting temperature AT = (7 — 1) T, 
where 7 is the bias factor and T is the temperature of the system. In the rest of the paper we will equivalently use 
either 7 or AT so as to simplify the notation. This parameter can be used to smoothly interpolate between unbiased 
sampling (7 = 1, AT = 0) and flat histogram (7 — 00 , AT —>• 00 ). One can thus introduce a set of replicas 
using different values of AT, ranging from 0 to a value large enough to allow all the relevant barriers to be crossed. 
Metadynamics relies on the accumulation of a history dependent potential and cannot be applied straightforwardly 
to a large number of CVs. In the next subsection we show that this issue can be circumvented by performing many, 
low-dimensional, concurrent metadynamics simulations. We then show how to combine many simulations of this kind 
in a multiple-replica scheme. 


Concurrent Well-tempered Metadynamics 

In well-tempered metadynamics a history dependent potential V (s, t) acting on the collective variable s is introduced 
and evolved according to the following equation of motion 

kaAT v(°.t) 

V{s,t) = -e K{s — s{t)) 

tb 

Here ks is the Boltzmann constant, T the temperature, tb is the characteristic time for the bias evolution, AT is a 
boosting temperature, and AT is a kernel function which is usually defined as a Gaussian. For simplicity we consider 
the case of a single CV. The variance of the Gaussian provides the binning in GV space and is usually chosen based 
on CV fluctuations or adjusted on the flv. |34) By assuming that the bias is growing uniformly with time one can show 
rigorously [ 2 II |35] that in the long time limit the bias potential tends to 

l^V{s,t) = -,^^F{s) + C{t) ( 1 ) 

so that the following probability distribution is sampled 

f(-) 

lim P{s,t) oc e 

t—>-oo 

The role of AT is that of setting the effective temperature for the CV. The explored conformations are thus taken 
from an ensemble where that CV only is kept at an artificially high temperature, similarly to other methods, |36H38j 
but has the nice feature that it is obtained with a bias that is quasi-static in the long lime limit. The bias is usually 
grown by adding a Gaussian every Nq steps. As a consequence, to obtain an initial growing rate equal to ^ 7 ^, the 
initial Gaussian height should be chosen equal to NqA t where At is the MD time step. 

We here propose to introduce a separate history-dependent potential on each CV 

fcoAT Vc(so,t) 

VM = — - e-^^K{s^ - s„(t)) 

Tb 
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where a = 1 ,... ,Ncv is the index of the CV and Ncv is the number of CVs. The growth of each of these bias 
potentials will depend only on the marginal probability for each CV 

P{Sa) oc 

(is\(is2 • ■ • dSa—1 ^ 2 , ■ ■ ■ 5 ^ Ncv ) 

In the long time limit, this potential will tend to flatten the marginal probabilities for every single CV. Since CVs are 
in general non orthogonal between each other, one should consider the fact that whenever a bias is added on a CV 
also the distribution of the other CVs is affected. In the following we will discuss this issue considering two CVs only, 
but the argument is straightforwardly generalized to a larger number of CVs. 

Two independent variables. If two CVs are independent, the joint probability is just the product of the two marginal 
probabilities, i.e. P(sq,,S/ 3 ) = Pa{sa)Pp{sp). Adding a bias potential on a CV will not affect the distribution of the 
other. As a consequence, in the long time limit the two bias potentials will converge independently to the predicted 
fraction of the free energy as in Eq. The final bias potential will be completely equivalent to that obtained 
from a two-dimensional well-tempered metadynamics, but will only need the accumulation of two one-dimensional 
histograms, thus requiring a fraction of the time to converge. A simple example on a model potential is shown in Fig. 
SI. 

Two identical variables. We also consider the case of two identical CVs, Sa = sp. This can be obtained for instance 
by biasing twice the same torsional angle. Here the potentials Va and Vp will grow identically, and the total bias 
potential acting on Sa will be Vtot = 214,. The total potential will grow as 


• 2fcijAr 

Vtot{Sa) = -e 


tb 

2kB^T vtot(-c,t) 


A"(Sq, ^a(t)) 


tb 


e X{Sa — Sa{t)) 


Thus, the net effect will be exactly equivalent to that of choosing a doubled AT parameter. In other words, the 
AT parameter acts in an additive way on the selected CVs. A similar effect can be expected if two CVs are linearly 
correlated. 

In realistic applications one can expect the behavior to be somewhere in the middle between these two limiting cases. 
The most important consideration here is that the bias potentials will tend to flatten all the marginal probabilities, 
but there will be no guarantee that the joint probability is flattened. Results for a simple functional form can be seen 
in Fig. S2. In the same figure it is possible to appreciate importance of using a self-consistent procedure when CVs 
are correlated. In ref |39j two metadynamics were applied on top of each other, namely on the potential energy and 
on selected CVs, in a non self-consistent way. This was possible because the correlation between the potential energy 
and the selected CVs is small. The need for a self-consistent solution was also pointed out in a recent paper [ID] where 
a generalization of the ABF method[22] was introduced. In that work independent one-dimensional adaptive forces 
were applied at the same time to different CVs so as to enhance the sampling of a high multidimensional space. 

In short, the novelty of the introduced procedure is that many low-dimensional metadynamics potentials are grown 
instead of a single multi-dimensional one. This allows the bias to converge very quickly to a flattening potential, 
with the degree of flatness controlled by the parameter AT. The flattening is expected to enhance conformational 
transitions which are otherwise hindered by free-energy barriers on the biased CVs. When variables are correlated 
the exact relationship between bias and free energy (Eq. could be lost. 


Hamiltonian Replica Exchange 

The procedure introduced above produces conformations in an ensemble which is in general difficult to predict. 
However, since the bias potential is known, one can in principle reweight results so as to extract conformations in 
the canonical ensemble. In the case of static bias potentials acting on the CVs, this can be done by weighting each 

frame as e ’‘b't . This can provide in principle correct results even if the joint probability is not flattened. It must 

be noticed that such a reweighting can provide statistically meaningful results only for small fluctuations of the total 
biasing potentials, on the order of fcsT. m However, in a typical setup one would be interested in biasing all the 
torsional angles of a molecule. Even if each of them contributes with a few fcsT, the total fluctuation of the bias 
would grow with the system size. For similar reasons, also the ABF-based scheme introduced in ref |40| is limited to 
a relatively low number of CVs. 
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A more robust and scalable procedure can be designed by introducing a ladder of replicas with increasing values 
of AT, ranging from 0 to a value large enough to enhance the relevant conformational transitions. The first replica 
(7 = 1, AT = 0) can be used to accumulate unbiased statistics. Replicas other than the first one feel multiple biasing 
potentials on all the CVs. From time to time an exchange of coordinates between neighboring replicas is proposed 
and accepted with probability chosen so as to enforce detailed balance with respect to the current biasing potential: 

acc = min (l, e~^) 

keT 

keT 

Here the sufhx i = 1,..., N^ep indicates the replica index, N^ep being the number of replicas. The exchanges allow 
the bias potential of every single replica to grow as close as possible to equilibrium taking advantage of the enhanced 
ergodicity of the more biased replicas. We notice that to reach a quasi-static distribution it is necessary that all 
the bias potentials converge for all the replicas. Since the time scale for convergence is related to the parameter 
tbjIH] it is convenient to use the same tb for all the replicas or, equivalently, to choose the initial deposition rate 
as proportional to AT. The number of replicas required to span a given range in the AT parameter is proportional 
to VAcv, thus allowing for a very large number of CVs. We underline that, if a very large number of CVs has to 
be calculated at every step on every replica, the overall performances could be degraded. This issue could be tackled 
using e.g. multiple time-step schemes. mm!] In the example presented here this performance issue was not observed. 

We notice that in principle one could use the bias potentials built with this protocol to perform a replica-exchange 
umbrella sampling simulation. In this manner the final production run would be performed with an equilibrium 
replica-exchange simulation. However, we observe that well-tempered metadynamics is designed so that the speed at 
which the bias grows decreases with time and the potential becomes quasi-static. In the practical cases we investigated, 
this second stage was not necessary. 


Model systems 


Alanine dipeptide 

Alanine dipeptide (dALA) was modeled with the Amber99SB-ILDN[43l|44j force field and solvated in an truncated 
octahedron box containing 599 TIP3P[1S] water molecules. The LINGS gaiiT] algorithm was used to constrain all 
bonds and equations of motion were integrated with a timestep of 2 fs. For each replica the system temperature was 
kept at 300 K by the stochastic velocity rescaling thermostat. [48] For all non-bonded interactions the direct space 
cutoff was set to 0.8 nm and the electrostatic long-range interactions were treated using the default particle-mesh 
EwaldjlH] settings. All the simulations were run using GROMACS 4.6.5(30] patched with the PLUMED 2.0 plugin. [ST] 
We underline that the possibility of running concurrent metadynamics within the same replica is a novelty introduced 
in PLUMED 2.0. 

The REGT simulation was performed with 6 replicas. The backbones dihedral angles (dt and 4)) and the gyration 
radius {Rg) were selected as GVs. The 7 factors were chosen from 1 to 15 following a geometric distribution. We 
recall that a geometric replica distribution is optimal for constant specific-heat systems. In REGT, this would be true 
if the exploration of each of the biased GVs were limited to a quasi-parabolic minimum in the free-energy landscape. 
Whereas this is clearly not true in real cases (e.g. double-well landscapes) we found that a geometric schedule was 
leading to a reasonable acceptance in the cases investigated here. The possibility of optimizing the replica ladder is 
left as a subject for further investigation. For the dihedral angles the Gaussian width was set to 0.35 rad and for the 
Rg to 0.007 nm. The Gaussians were deposited every 500 steps. The initial Gaussian height was adjusted to the AT 
of each replica, according to the relation h = Nc^t, in order to maintain the same tb = 12 ps across the entire 
replica ladder. The GVs were monitored every 100 steps, and exchanges were attempted with the same frequency. 
The simulation was run for 20 ns per replica. 

A H-REMD simulations where the force-field dihedral terms were scaled (H^i/i-REMD) was also performed, as 
implemented in an in-house version of the GROMAGS code. [52] The same initial structures, number of replicas and 
simulated time as in REGT was used. The scaling factor A for each replica was selected using the relation A = I /7 to 
allow for a fair comparison of REGT and H-REMD. Finally a conventional MD simulation in the NVT ensemble was 
run for 120 ns using the same settings. 
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FIG. 1: Schematic representation of the collective variables used for the tetranucleotide simulation. For each nucleotide, the 
labeled dihedral angles and the minimum distance between each nucleobase center of mass and the other three nucleobases 
were biased. 


Tetranucleotide 

The second system considered was an RNA oligonucleotide, sequence GACC. The initial coordinates were taken from 
a ribosome crystal structure (PDB: 3G6E), residue 2623 to 2626. Simulations were performed using the Amber99- 
bsc 0 ,^OL 3 force field. [32 ESI EH The system was solvated in a box containing 2502 TIP3P|3S] water molecules and the 
system charge was neutralized by adding 3 Na+ counterions, consistently with previous simulations. [Ml ES] A REGT 
simulation was performed using 16 replicas simulated for 300 ns each. The 7 ladder was chosen in the range from 1 
to 4 following again a geometric distribution. The initial structures for the H-REMD were taken from a 500 ps MD 
at 600 K, to avoid correlations of the bias during the initial deposition stage of the WT-MetaD. Other details of the 
simulation protocol were chosen as for the previous system. As depicted in Fig. [l] for each residue the dihedrals of 
the nucleic acid backbone (a, /3, e, 7, ir ), together with the pseudo-dihedrals angles of the ribose ring ( 0 i and 62 ) 
and the glycosidic torsion angle (x) were chosen as GVs. To help the free rotation of the nucleotide heterocyclic base 
around the glycosidic bond, the minimum distance between the center of mass of each base with the other three bases 
was also biased. For the WT-MetaD we used the same parameters as in the previous system. Gaussian width for the 
minimum distance between bases was chosen equal to 0.05 nm. 

For this system a H^i/t-REMD, a T-REMD and a plain MD simulation were performed in addition to the REGT. In 
the case of H^i/i-REMD we used 24 replicas with scaling factors A ranging from 1 to 0.25, so as to cover the same range 
of the 7 values chosen for the REGT. In the T-REMD 24 replicas were used to cover a temperature range between 300 
K and 400 K with a geometric distribution. For both methods, T-REMD and H^i/j-REMD, the simulation length was 
200 ns per replica. Exchanges were attempted every 120 steps. The conventional MD simulation was run for 4.8 /rs. 
All the simulations (REGT, H^i/j-REMD, T-REMD, and conventional MD) correspond to the same total simulated 
time. 


Analysis 


Dihedral entropy 

As the bias compensates the underlying free energy the probability distribution of the biased GVs is partially 
flattened. The main GVs used in our method are dihedrals angles. To quantify the effect of the Hamiltonian 
modifications on the angle distributions one-dimensional entropies {Sid) were estimated. The calculation procedure 
was equivalent to the one used in ref (SHI to evaluate the configurational entropy associated with soft degrees of freedom 
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in proteins. We employed wrapped Gaussian kernels to estimate the histogram profile of each dihedral. Histograms 
were calculated with PLUMED 2.0. For all the distributions the bandwidth for the kernel density estimation was set 
to 0.017 rad. We underline that using this definition we only evaluate the flatness of the individual one-dimensional 
distributions, and cross-correlation between CVs is ignored. 


RNA conformations 

RNA conformations were classified according to the combination of the y angle rotameric state of each nucleotide. 
Torsion orientations in the range of -0.26 to 2.01 rad were consider as syn, while the remaining ones were classified 
as anti. The limiting values were chosen according to the position of the barriers in the y free-energy profiles of all 
the residues. The result of this clustering procedure gave 2'^ = 16 different states that are kinetically well separated 
by the high torsional barriers. We observe that the population of these states does not depend only on the torsional 
potential associated to the y dihedrals but include contributions from base-base stacking, hydrogen bonds, solvation 
of bases, etc. 


Results 

In this section we first show the validation of our methodology on a standard model system, dALA in water. 
Then we present results for the more challenging case of the conformational sampling of a tetranucleotide. For all 
the applications we benchmark against plain MD and a H-REMD where the dihedral potentials are scaled. All the 
comparisons are made using the same total simulated time. 


Alanine Dipeptide 

The goal of the introduced method is to enhance conformational sampling in the unbiased replica. The possibility 
to explore different metastable conformations in this replica relies on the fact that probability distributions in the 
biased replicas are flattened and that conformations can travel across the replica ladder. These conditions can be 
verified by monitoring the exchange rate and the flatness of the distributions. 

The acceptance rate is in the range 65-72% for RECT and in the range 43-53% for Hj^^^i-REMD, indicating that 
the former method requires less replicas. This is likely due to the fact that the total number of scaled dihedrals in 
Hdi/i-REMD is larger than the number of biased CVs in RECT. For both REMD methods we also verified that all 
the trajectories in the generalized ensemble sampled the same conformational ensemble (see Fig. S3). 

A quantitative measure of the flatness of the distribution in the biased replicas can be obtained from the dihe¬ 
dral entropy, shown in Fig. as a function of scaling factors (7 and A for RECT and H^i/t-REMD, respectively). 
The limiting value corresponding to a flat distribution is also indicated. Entropy grows faster as a function of the 
scaling factor when using RECT, indicating that free-energy barriers on the dALA isomerization transition are more 
effectively compensated by the bias potentials. With Hj/i^^-REMD entropy of 4/ angle saturates and apparently the 
distribution cannot be further flattened by decreasing A. In the case of the $ angle, the dihedral entropy does not 
grow monotonically when A is decreased. This behavior indicates that the relevant free-energy barriers are not only 
originating from the dihedral force-field terms. The conformational transitions involve indeed also changes in water 
coordination, reorganization of hydrogen bonds, non-bonded interactions, etc. On the contrary, RECT achieves an 
almost flat distribution for both dihedral angles at the highest value of the 7 factor. Backbone dihedral distributions 
for all the replicas are shown in Fig. S4. The conformations sampled on each replica are shown projected on the $,'1' 
free-energy landscape in Fig. S5, where it can be appreciated that all the relevant basins (a, /3, and a^j) are explored 
and connected by points close to the minimum-action pathways, (see refs gniisT]). 

To assess the efficiency and the accuracy of the introduced enhanced sampling technique the free-energy difference 
AF between the states cj) € [—tt, 0] and (f) S [0, |] was calculated from the distribution of the unbiased replica. Results 
are shown as a function of time in Fig. for the two REMD schemes and for the reference conventional MD. Both 
H-REMD methods converge to the right value with a similar behavior, whereas plain MD needs several tens of ns for 
the first transition to be observed. The similarity in the convergence of RECT and H^i/j-REMD indicates that for 
this system the moderate flattening of the distribution induced by H^ih-REMD is sufficient to achieve ergodicity on 
this time scale. In order to better evaluate differences between the performance of RECT and H^^i^-REMD we applied 
this methodologies to a more complex system. Results are shown in the next section. 
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FIG. 2: Entropy for (top) and (bottom) dihedral angles in alanine dipeptide. Entropies are shown as a function of 1/A 
and 7 for H^ih-REMD and RECT respectively. As the entropy increases the dihedral distributions become more flat. The 
maximum entropy value corresponding to a flat distribution is represented with a straight line. 



FIG. 3: Estimate of the free-energy difference between the two metastable minima in alanine dipeptide. Data are shown for 
both replica-exchange methods (Hdi/i-REMD and RECT) and for conventional MD as a function of the total simulated time. 


Tetranucleotide 

Also in this case we monitor the average exchange ratio (76-83% for H£iift,-REMD, 25-32% for T-REMD, and 60-80% 
for RECT). In Fig. S7 the variation of the exchange ratios in time is shown for the exchanges between the first 2 
and the last 2 replicas of each method. We also checked the consistency of trajectories along the replica ladder. As 
it can be appreciated in Fig. S6, for Hdift,-REMD and RECT the empirical distribution of RMSD is very similar 
for all the trajectories in the generalized ensemble, indicating that, for each method, all of them sampled the same 








conformational space. On the contrary, in the case of T-REMD, agreement among the distributions of RMSD is 
very poor. During this simulation trajectories across the temperature space remain trapped on different metastable 
conformations. The same behavior was obtained in ref [5^ were several T-REMD simulations were performed on the 
same system, with the same number of replicas and a similar temperature range. In that work divergence among the 
obtained generalized ensembles was observed even for a simulated time as long as 2 fis per replica. For H^i/i-REMD 
and RECT round-trip times are shown on Fig. S 8 . The average round-trip time is « 0.5 ns for H^i/j-REMD, « 1.8 
ns for T-REMD, and « 1.2 ns for RECT. 

In Fig. I^we show the sum of the entropies for the 32 dihedrals used as CVs. In this respect, RECT is clearly more 
effective than H^i/i-REMD in flattening the dihedral distributions, consistently with what was observed for dALA. 
Notably, the entropic increment observed in RECT is close to the one observed in T-REMD when using an equivalent 
temperature. This confirms that RECT has an effect comparable to that of raising the temperature of the biased 
CVs by a factor 7 . 

The significance of this entropic values could be appreciated on the time series and related histograms for all the 
dihedral angles shown in Fig. S9-12 for the most and least ergodic replica of H^i^t-REMD and RECT. It is clear that 
for RECT, at the most ergodic replica, all the accessible torsional range is sampled. On the contrary, in the highest 
replica of H^i/i-REMD the distributions of some torsions are not flattened. 


l/A 



FIG. 4: Total entropy of backbone, puckering and glycosidic dihedral angles in the tetranucleotide for both replica-exchange 
methods. Entropies are shown as a function of 1/A and 7 , for Hdi/i-REMD, RECT respectively. Eor T-REMD, temperature 
is chosen as T = ySOOA. As the entropy increases the dihedral distributions become more flat. The maximum entropy value 
corresponding to a flat distribution is represented with a straight line. Entropies obtained for the unbiased replicas in the three 
methods are consistent within their error bars (error not shown). 


The transition around the glycosidic bond, from anti to sj/n, is among the slowest relaxation times in RNA 
dynamics. |58| To evaluate the convergence of the unbiased replica we analyzed the population of the anti rotamer 
for each nucleotide x angle. Populations are shown in Fig. as a function of the total simulated time. For all the 
nucleotides the anti conformations are preferred. The guanosine is the nucleotide with the highest syn proportion, 
and the cytidines the ones with the smallest (< 2%), as correspond to their rotameric preferences. [55] Values from 
both H-REMD approaches seem well consistent, except for the population of the first nucleotide. From the time 
behavior of these populations, it is clear that for all the REMD approaches the guanosine proportion of anti is the 
most difficult to converge. Here RECT can reach values close to a longer reservoir-REMD simulation|5^ while both 
Hdi/i-REMD and T-REMD show results closer to those obtained from conventional MD, with a higher occupation of 
the anti conformer. 


We observe that our method is enforcing the exploration of both anti and syn conformations in the biased replicas 
for each nucleotide independently. This however does not guarantee that all the 16 combinations of anti and syn 
conformations are explored. Fig. shows the free energy of the RNA structures grouped by the combination of the 
X angle anti(a)/syn(s) rotamers. All 16 combinations, except for ssss and asss, are sampled in the unbiased replica 
from RECT. On the contrary, the unbiased replica from T-REMD and H^ih-REMD explores respectively 13 and 8 
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FIG. 5; Estimated glycosidic angle anti population for each nucleotide as a function of the total simulated time. Data are 
shown for H^iji-REMD, T-REMD, and RECT unbiased replicas and for conventional MD. Reference values taken from ref |26| 
are shown as dashed lines. 


of the states, and plain MD only 5 of them. The most populated cluster corresponds to an all-anti conformation, 
followed by the saaa. Then, the three clusters asaa, ssaa, and sasa appear with similar population. 

In the same figure the free-energy values for the ergodic replica show that all the 16 combinations are populated 
in RECT within a range of QksT. In the case of H^i/j-REMD the most ergodic replica visits only 9 combinations 
with a population that is very close to that of the unbiased replica. The most ergodic replica in T-REMD explores 14 
clusters, but their populations have a large statistical errors. We highlight the fact that results from T-REMD could 
be affected by the lack of convergence of trajectories across the temperature space (see Fig. S 6 ). This could lead to 
an underestimation of the errors as evaluated from block analysis. 


Discussion 

The introduced method allows to build bias potentials for a Hamiltonian replica-exchange scheme using concurrent 
well-tempered metadynamics on several CVs. Replicas are simulated using a ladder of well-tempered bias factors 7 . 
When CVs are correlated, the self-consistency among the bias potentials is crucial to achieve flat sampling in each 
individual CV, as illustrated in Fig. S2. In this case the exact relationship between bias and free energy is lost. 
We also remark that here flattening is not complete but modulated by the value of 7 . This is useful since it avoids 
sampling very high energy states (e.g. with steric clashes) that would have a very low chance of being accepted in 
the unbiased replica. The method compares favorably with both conventional MD and H^i/j-REMD. The method 
slightly outperforms T-REMD, where the entire system is heated, indicating that for these small systems there is not 
a substantial advantage in schemes where part of the system is biased. However, RECT can be straightforwardly 
generalized to large systems since the acceptance only depends on the size of the biased portion. 

Results from both dALA and tetranucleotide simulations show that the bias potentials constructed with concurrent 
WT-MetaD are able to gradually scale the free-energy barriers. We notice that only barriers in the one-dimensional 
free-energy profiles are compensated, which means that some regions in the multidimensional space of all the CVs 
might not be explored. In principle this could hide some important minima that would never be observed. We did 
not observe this problem in the applications presented here. 

The second application on which we tested RECT, namely conformational sampling of a tetranucleotide, is par¬ 
ticularly challenging. The conformational space of these small RNA molecules is not constrained by Watson-Crick 
pairings and ergodic sampling is out of reach of conventional MD simulations. |26L 155] So far, converged ensembles have 
been obtained only trough highly expensive multidimensional REMD simulations, corresponding to a total simulated 
time of several tens of fis.|181 m\ One of the reasons for this difficult convergence is the long relaxation time for 
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FIG. 6: Estimated free energies for the tetranucleotide conformations clustered according to the x angle anti/syn rotameric 
combinations (circles). Free energies are computed as —ksTlogPi, where Pi is the normalized population of each cluster on 
the indicated replica. Grey boxes represent relative populations higher than 1%. Gonfidence intervals are shown as bars and 
span the range [—kBT\og{Pi + APi), —fcsTlog(Pi — APi)], where APi is the standard deviation of the average Pi as obtained 
from four blocks. Clusters which are observed in only one of the four blocks have an infinite upper bound. 


the anti to syn transitions, which could be additionally hindered by an incorrect force-field description of base-base 
stacking and base-solvent interactions. [55] 

Fig. [^illustrates the ability of RECT to accelerate conformational transitions among the y angle anti/syn rotamers. 
Although the conformational space of the more biased replicas is highly expanded, the convergence in the unbiased 
replica is not affected. On the contrary the method facilitates the sampling of glycosidic rotamer conformations that 
otherwise would not be explored by MD simulations of the same overall length. We finally remark that our procedure 
can be combined with weighted histogram[50] so as to include the statistics of the biased replicas. 


Comparison with related state-of-the-art methods 

RECT is based on the idea of building a replica ladder where a large set of selected CVs is progressively heated. 
CVs are heated by flattening their distribution with concurrent well-tempered metadynamics. We first discuss the 
possibility of using methods other than well-tempered metadynamics to build the replica ladder. Possible alternatives 
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here include ABF [22] or a recently proposed variational approach inn- These methods could be used in a RECT 
scheme provided they are suitably extended so as to sample a partially flattened distribution. We also observe 
that other methods aimed at keeping selected CVs at a given temperature have been proposed based on coupling 
thermostats to CVs directly. [5BH55] These techniques have been mostly used in the past with an exploration purpose 
relying on additional calculations so as to provide free energies (see, e.g., ref. [H 2 ]) but it is not clear if they can be 
integrated in a RECT scheme. 

In the following we discuss the comparison of RECT with related methods that are not directly based on CV 
tempering. 

Comparison with H-REMD of Curuksu and Zacharias. Our method is closely related to the one introduced in ref 
m- There, a bias potential aimed at disfavoring the most probable rotamers is manually constructed and applied on 
several replicas using a scaling factor. This bias disfavors the major minima but does not ensure a proper compensation 
of the free-energy barriers, as their positions and magnitudes are not a priori known. The main advantage of RECT 
is that several low-dimensional bias potentials are built with a self-consistent procedure so that the technique can be 
straightforwardly applied to a large number of degrees of freedom. 

Comparison with bias-exchange metadynamics. In bias-exchange metadynamics every replica performs an indepen¬ 
dent metadynamics simulation so that one CV at a time is feeling the flattening potential. Thus, it is typically used 
with a relatively small number of ad hoc designed CVs capable to describe the relevant conformational transitions. On 
the other hand, RECT is designed to be used with a very large number of dummy CVs with little a priori information 
and to bias them concurrently to exploit their cooperation in enhancing conformational sampling. For this reason, 
the two approaches are complementary and could even be combined in a multidimensional replica exchange suitable 
for a massively parallel environment. 

Comparison with solute tempering and related methods. In replica exchange with solute tempering the solute 
Hamiltonian is scaled so as to obtain an effect equivalent to a rise in the simulation temperature. mm Any set of 
atoms can be identified as solute, giving the opportunity to enhance sampling in a region localized in space. [5^ 163] 
This requires modifying charges of the enhanced region, with long range effects and sometime affecting fundamental 
properties such as hydrophobicity. In our method, the bias potentials act on precisely selected degrees of freedom 
minimally perturbing their coupling with the rest of the system. Moreover, the bias is adaptively built so as to 
compensate the free energy and not the potential energy, so that with properly chosen CVs it could be used to 
compensate entropic barriers. 

Comparison with hyperdynamics and accelerated MD. In these methods the potential energy of the system is 
modified so as to decrease the probability to sample minima on the potential energy. On the contrary, 

RECT employs a bias which is related with the free energy so as to achieve a flatter histogram on the selected CVs. 

We finally remark that RECT, although formally based on the a priori choice of a set of CVs, typically requires 
the same amount of information of methods not based on CVs. Indeed, as we have shown, the method can be 
easily applied to a very large number of CVs, virtually including by construction all the slow degrees of freedom 
of the system. Additionally, when a few relevant CVs can be identified based on chemical intuition, RECT can be 
straightforwardly combined with standard metadynamics similarly to parallel tempering [ 66 ] or solute tempering. m 


Conclusion 

Replica exchange with collective-variable tempering (RECT) has been here proposed as a novel and flexible 
enhanced-sampling method. RECT takes advantage of the adaptive nature of well-tempered metadynamics to build 
bias potentials that compensate free-energy barriers. The flattening of the barriers is modulated by the well-tempered 
factor 7 , and the chosen collective variables (CVs) are effectively kept at a higher temperature. The biasing potentials 
are built combining concurrent low-dimensional metadynamics protocols so as to be usable on a very large number of 
CVs. Multiple replicas are then used so as to smoothly interpolate between a highly biased, ergodic simulation and 
an unbiased one (7 = 1). The number of required replicas scales with the square root of the number of chosen CVs 
for a fixed range of 7 factors. This allows a very large number of CVs to be biased, so that virtually all the relevant 
transitions can be accelerated. The CVs used here were mostly dihedral angles, which exhibit relevant barriers in 
many biomolecular conformational transitions, but the method can be used with any CV. The application of this 
technique to the dALA in water shows that the CV probability distributions are effectively flattened by the action 
of the bias potentials and unbiased statistics is correctly recovered. In the case of the tetranucleotide conformational 
sampling is greatly enhanced since RECT effectively overcomes the high free-energy barriers of the x angle transitions 
that hindered the conformational sampling at room temperature. RECT is available in PLUMED and a sample input 
file is provided in Fig. S13. RECT is a promising tool to enhance the exploration of the conformational space in 
highly flexible biomolecular systems such as RNA, proteins, or RNA/protein complexes. 
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II. SUPPLEMENTARY INFORMATION 

Examples of biasing concurrently different CVs on model potentials with independent (Fig. SI) and correlated 
variables (Fig. S2). RMSD distributions of dALA from trajectories across the replica ladder (Fig. S3). Histograms 
of dALA dt and $ dihedral angles for all replicas (Fig. S4). The projection of each replica trajectory on the dihedral 
free-energy landscape (Fig. S5). RMSD distributions of RNA tetranucleotide from trajectories across the replica 
ladder (Fig. S6). Average exchange ratios vs time (Fig. S7). Round-trip times for the tetranucleotide simulations 
(Fig. S8). Time series and histograms of the torsional angles for the REMD lower and higher replicas (Fig. S9-S12). 
Example of a PLUMED 2.1 input file for a RECT simulation (Fig. S13). 
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